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Figure 7 
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Figure 8A 
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Figure 8B 
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Figure 9 
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Figure 10A 
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Figure 10B 
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Figure 11B 
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Figure 12 
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Figure 13B 
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SPLIT-REMERGE METHOD FOR 
ELIMINATING PROCESSING WINDOW 
ARTIFACTS IN RECURSIVE 
HIERARCHICAL SEGMENTATION 

5 

CROSS-REFERENCE TO RELATED 
APPLICATIONS 

The present application is a continuation in part of U.S. Ser. 
No. 10/845,419 filed May 11, 2004 now abandoned, and to 
which is incorporated herein by reference. 

ORIGIN OF THE INVENTION 

The invention described herein was made by an employee 15 
of the United States Government, and may be manufactured 
and used by and for the Government or for governmental 
purposes without the payment of any royalties thereon or 
therefor. 

20 

FIELD OF THE INVENTION 

The present invention relates to a method for eliminating 
processing window artifacts that may occur when partition- 
ing spatially related data into sections or regions for process- 25 
ing. More particularly, the present invention describes a split- 
remerge method for identifying and splitting possible 
artifacts of data and then merging that data with the most 
appropriate region. 

30 

BACKGROUND OF THE INVENTION 

Segmentation, the partitioning of data into related sections 
or regions, is a key first step in a number of approaches to data 
analysis and compression. In data analysis, the group of data 35 
points contained in each region provides a statistical sampling 
of data values for more reliable labeling based on data feature 
values. In data compression, the regions form a basis for 
compact representation of the data. The quality of the prereq- 
uisite data segmentation is a key factor in determining the 40 
level of performance of most of these data analysis and com- 
pression approaches. 

Most segmentation approaches can be placed in one of 
three categories: 

(i) characteristic feature thresholding or clustering, 45 

(ii) boundary detection, or 

(iii) region growing. 

Characteristic feature thresholding or clustering does not 
exploit spatial information, and thus ignores information that 50 
could be used to enhance the segmentation results. While 
boundary detection does exploit spatial information by exam- 
ining local edges found throughout the data, it does not nec- 
essarily produce closed connected region boundaries. For 
simple noise-free data, detection of edges usually results in 55 
straightforward region boundary delineation. However, edge 
detection on noisy, complex data often produces missing 
edges and extra edges that cause the detected boundaries to 
not necessarily form a set of close connected curves that 
surround connected regions. Region growing approaches to 60 
segmentation are preferred because region growing exploits 
spatial information and guarantees the formation of closed, 
connected regions. 

Segmentation is often used in the analysis of imagery data. 
The techniques described can be applied to image data and to 65 
any other data that has spatial characteristics. A data set has 
spatial characteristics if it can be represented on an n-dimen- 
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sional grid, and when so represented, data points that are 
nearer to each other in the grid generally have a higher sta- 
tistical correlation to each other than data points further away. 
For remotely sensed images of the earth, an example of a 
segmentation would be a labeled map that divides the image 
into areas covered by distinct earth surface covers such as 
water, snow, types of natural vegetation, types of rock forma- 
tions, types of agricultural crops and types of other man 
created development. In unsupervised image segmentation, 
the labeled map may consist of generic labels such as region 
1 , region 2, etc., whichmay be converted to meaningful labels 
by a post- segmentation analysis. In image analysis, the group 
of image points contained in each region provides a good 
statistical sampling of image values for more reliable labeling 
based on region mean feature values. In addition, the region 
shape or texture can be analyzed for additional clues to the 
appropriate labeling of the region. 

A segmentation hierarchy is a set of several segmentations 
of the same data at different levels of detail in which the 
segmentations at coarser levels of detail can be produced 
from simple merges of regions at finer levels of detail. This is 
useful for applications that require different levels of segmen- 
tation detail depending on the particular data objects seg- 
mented. A unique feature of a segmentation hierarchy that 
distinguishes it from most other multilevel representations is 
that the segment or region boundaries are maintained at the 
finest data granularity for all levels of the segmentation hier- 
archy. 

In a segmentation hierarchy, an object of interest may be 
represented by multiple segments in finer levels of detail in 
the segmentation hierarchy, and may be merged into an 
encompassing region at coarser levels of detail in the segmen- 
tation hierarchy. If the segmentation hierarchy has sufficient 
resolution, the object of interest will be represented as a single 
region segment at some intermediate level of segmentation 
detail. The segmentation hierarchy may be analyzed to iden- 
tify the hierarchical level at which the object of interest is 
represented by a single region segment. The object may then 
be identified through its spectral and region characteristics, 
such as shape and texture. Additional clues for object identi- 
fication may be obtained from the behavior of the segmenta- 
tions at the hierarchical segmentation levels above and below 
the level at which the object of interest is represented by a 
single region. 

In U.S. Pat. No. 6,895,1 15, which is incorporated herein by 
reference, a segmentation approach is described that auto- 
matically provides hierarchical segmentations for data at sev- 
eral levels of detail. This approach, called HSEG, is a hybrid 
of region growing and spectral clustering that produces a 
hierarchical set of segmentations based on detected natural 
convergence points. Because of the inclusion of spectral clus- 
tering, the HSEG algorithm is very computationally inten- 
sive, and cannot be performed in less than a day on moder- 
ately sized data sets, even with the most powerful single 
processor computer currently available. The processing time 
problem was addressed through a recursive formulation of 
HSEG, called RHSEG. RHSEG can process moderately 
sized data sets in a reasonable amount of time on currently 
available PCs and workstations. Larger data sets required the 
use of a parallel implementation of RHSEG on a parallel 
computing system. 

However, a problem with the RHSEG algorithm and cer- 
tain other data processing algorithms that similarly subdivide 
and subsequently recombine data during processing instant 
processing artifacts can be introduced by the division and 
recombination of the data. An example of these processing 
artifacts can be demonstrated on an 896x896 pixel section of 
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Landsat ETM+ (Enhanced Thematic Mapper) data displayed 
in FIG. 1. This image was obtained on May 28, 1 999 over the 
southwestern coast of the eastern shore of Maryland. The six 
non-thermal bands were used in the segmentation tests. 

FIG. 2 displays the segmentation of this image into 96 5 
regions as produced by the basic version of RHSEG. Straight- 
line processing window artifacts from the recursive quarter- 
ing of the image data are quite evident in sub-regions of the 
Chesapeake Bay 1,2 and on land in the right center part of the 
image 3 . These straight or blocked lines have no relationship to 
to the data in the image itself, but arose from the recursive 
subdivision and recombination process. Processing artifacts 
may be very noticeable as lines, or merely a few mislabeled 
pixels. 

In the prior art, a “contagious clusters” or “contagious 15 
regions” concept has been used to attempt to reduce or elimi- 
nate the processing window artifacts. The contagious regions 
concept can be described as follows: 

Flag any region that touches a boundary between process- 
ing windows and suppress any merging between flagged 20 
regions and any other region. 

If a non-flagged or “non-contagious” region attempts to 
merge with a flagged or “contagious” region, the previously 
non-flagged region becomes flagged or “contagious.” 

Thus, the contagious property of the flagged regions is 25 
literally contagious. Unfortunately, when the contagious 
regions concept is applied to the RHSEG algorithm, and 
when more than two or three levels of recursion are utilized, 
the RHSEG algorithm is only able to effectively process the 
data. The RHSEG algorithm effectively stalls because so 30 
many regions become contagious that the number of regions 
in the processing window becomes so large that the process- 
ing time required is not sufficiently advantageous over a 
non-segmented image processing approach. 

Indirect mechanisms, however, also exist. For example, 35 
increasing the number of regions at which convergence is 
achieved at intermediate levels of the recursive processing 
may indirectly cause a reduction in processing window arti- 
facts. A larger value may delay some region merging deci- 
sions that would have involved regions on the borders of 40 
processing windows to occur after those regions are no longer 
on the borders of processing windows. This indirect method, 
however, is inefficient because processing time increases with 
larger values of the number of regions needed to achieve 
convergence. Further, processing artifacts are not always 45 
eliminated via this method. Other approaches to reducing 
window artifacts may manipulate other parameters in the 
recursive hierarchical segmentation processing, but also 
increase processing time and resources, such that the 
approaches become impractical for large sets of data. 50 

Indeed, all previously developed techniques for splitting 
inappropriately merged pixels or processing image data in a 
fashion to avoid creating window artifacts unacceptably 
increase the processing time required. Thus, in prior applica- 
tion Ser. No. 1 0/845,4 1 9, a switch-pixels method of address- 55 
ing window artifacts was disclosed, however, this technique 
had no mechanism for giving priority to spatial adjacency in 
switching pixels from one region to another. 

SUMMARY OF THE INVENTION 60 

Accordingly, it is an object of the present invention to 
implement a split-remerge process for eliminating processing 
window artifacts in recursive hierarchical segmentation of 
data. The foregoing object of the invention is achieved by 65 
identifying candidate pixels or data points or regions of points 
that have been inappropriately merged, identifying a candi- 
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date region with which those data points or sets might be more 
appropriately merged and then evaluating the best merger 
candidate giving due weight to spatial distance between the 
flagged data and the candidate region for merger. While ini- 
tially designed for the analysis of single-band, multi spectral 
or hyperspectral remotely sensed imagery data for earth sci- 
ence applications, the software innovation also has applica- 
tions to image data compression of image data archives, data 
mining (searching for particular shapes of objects with cer- 
tain feature vector characteristics), and data ffision (based on 
matching region features between data sets from different 
times and/or different sensors). Applications outside of 
remote sensing are the analysis of imagery for medical appli- 
cations and for nondestructive evaluation in manufacturing 
quality control. A possible military application is land mine 
detection. 

BRIEF DESCRIPTION OF THE DRAWINGS 

These and other advantages of the invention will be more 
fully appreciated from the following description of preferred 
embodiments taken in conjunction with the accompanying 
drawings, of which: 

FIG. 1 is an example of a satellite image before segmenta- 
tion; 

FIG. 2 is an example of the satellite image of FIG. 1 after 
segmentation into 96 regions according to the basic RHSEG 
algorithm with spclust_wght=l .0. 

FIG. 3 is an example of the satellite image of FIG. 1 after 
segmentation into 8 regions and utilizing the switch pixels 
approach to eliminating processing artifacts and 

spclust_wght=0. 1 . 

FIG. 4 is an example of the satellite image of FIG. 1 after 
segmentation into 8 regions and utilizing the switch pixels 
approach to eliminating processing artifacts and 

spclust_wght=l .0. 

FIG. 5 is an example of the satellite image of FIG. 1 after 
segmentation into 10 regions by the RHSEG algorithm with 
the spclust_wght=l .0 and utilizing the split-remerge method 
of eliminating processing artifacts according to the invention. 

FIG. 6 is an example of the satellite image of FIG. 1 after 
segmentation into 3 1 regions using the split -remerge proce- 
dure of the invention to process artifacts, but without the 
contagious pixel aspect implemented. 

FIG. 7 is a flow chart illustrating the HSWO algorithm. 

FIGS. 8A and 8B are a flow chart illustrating the HSEG 
algorithm. 

FIG. 9 is a flow chart representing a recursive formulation 
of the HSEG algorithm. 

FIGS. 10A and 10B are a flow chart representing the 
“switch pixels” rhseg(level,X) algorithm. 

FIGS. 11A and 1 IB are a flow chart representing the “split- 
remerge” rhseg(level,X) algorithm. 

FIG. 12 is a flow chart representing the recursive remerge 
(level, max_threshold, X) function. 

FIGS. 13A and 13B are a flow chart representing the 
restricted HSEG algorithm. 

FIG. 14 is a graphical representation of the relationships 
between the parallel tasks executing the RHSEG program on 
a parallel computer for a 2-spatial dimension input data set. 

DETAILED DESCRIPTION OF THE INVENTION 

Most prior techniques for region growing data segmenta- 
tion are based on the classic definition of 2-spatial dimension 
image segmentation: 
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Let X be a 2-dimensional array representing an image. A 
segmentation of X can be defined as a partition of X into 
disjoint subsets X 1? X 2 , . . . , X^, such that 


1) U Xi = x 

i = 1 

2) Xi, i = 1, 2, ... , N is connected. 

3) P{Xi ) = TRUE for / = 1, 2, ... ,N, and 10 

4) P(Xi U Xj) = FALSE for i * j, 
where X t and Xj are adjacent 

15 

P(X Z ) is a logical predicate that assigns the value TRUE or 
FALSE to X z , depending on the image data values in X z . 

For example, a logical predicate based on the vector 
2 -norm is: on 


PiXi ) = 


^ X \\xj-xi\\ 2 z r ‘ 


( 1 ) 

25 


where n. is the number of pixels in region i, x.is the mean 
vector for region i, and T is a threshold. 

30 

These conditions might be summarized as follows: The first 
condition requires that every picture element (pixel) must be 
in a region. The second condition requires that each region 
must be connected, i.e. composed of contiguous image pixels. 
The third condition determines what kind of properties each 35 
region must satisfy, i.e. what properties the image pixels must 
satisfy to be considered similar enough to be in the same 
region. The fourth condition specifies that, in the final seg- 
mentation result, any merging of any adjacent regions would 
violate the third condition. If classical image segmentation is 40 
not implemented recursively, it will not be subject to prob- 
lems with processing window artifacts. A window is one of 
the data subsections defined by the recursive subdivision and 
subsequent recombination of the data in the context of “a 
processing window.” Processing Window Artifacts are sub- 45 
optimal region assignments for data elements caused by the 
structure of the recursive subdivision and subsequent recom- 
bination of the data in the context of a recursive implemen- 
tation of a segmentation algorithm. A recursive implementa- 
tion of a segmentation algorithm can also be referred to as a 50 
recursive segmentation algorithm. A recursive implementa- 
tion of a segmentation algorithm is defined as follows: (a) 
recursively divide the data into subsections until each subsec- 
tion consists of no more than a predetermined maximum 
number of data elements; (b) perform the segmentation algo- 55 
rithm on each data subsection such that no more than a pre- 
determined number of regions is reached. Then return the 
segmentation result and maximum merging threshold to the 
calling recursive level, (c) After a return from a deeper level of 
recursion, initialize the segmentation algorithm with segmen- 60 
tation results from all the deeper levels of recursion, and 
perform the segmentation algorithm on the data subsection at 
the current level of recursion such that no more than a prede- 
termined number of regions is reached, (d) After step (b) is 
completed for all data subsections from step (a), and after step 65 
(c) is completed for all recursive returns, continue the seg- 
mentation algorithm to completion on the entire data set. 
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A problem with the classic definition of image segmenta- 
tion is that the segmentation so defined is not unique. The 
number, N, and shape of the partitions, X 1; X 2 , . . . , X^, 
depend on the order in which the image pixels are processed. 
In addition, there is no concept of optimality contained in this 
definition of image segmentation. Under this classic defini- 
tion, all partitions that satisfy the conditions represent equally 
good or valid segmentations of the image. 

A less commonly used, more computationally intensive 
but often more accurate approach to 2-spatial dimension 
image segmentation is the Hierarchical Stepwise Optimiza- 
tion (HWSO) algorithm of Beaulieu and Goldberg. HSWO 
can be defined recursively as follows: 

Let X be a 2-dimensional array representing an image and 
let X l5 X 2 , . . . , X^ , X^be a partition of X into N regions 
such that 


N 

1) U Xi = X and 
/= 1 

2) Xi, i = 1, 2, ... , N is connected. 


Let G(X Z ) be a function that assigns a cost to partition X z , 
depending on the image data values in X z . Reorder the 
partition X 19 X 2 , . . . , X^, X^ such that G(X JV _ 1 UX JV )i 
G(X z UX 7 ) for all fyj where X^ and X^are adjacent and 
X z - and Xj are adjacent. The segmentation of X into N-l 
regions is defined as the partition X\, X' 2 , . . . , Xfy^ 
where X',=X, for i=l, 2, . . . , N-2 andXfy.^X^UX^. 

The initial partition may assign each image pixel to a separate 
region, in which case the initial value of N is the number of 
pixels in the image (N p ). Any other initial partition may be 
used, such as a partitioning of the image into nxn blocks, 
where n 2 «N^. 

For example, a cost function based on the vector 2-norm is: 


CTO 


- y \\xj-xi \\ 2 , 

n ‘ x ;<= Xj 


(2) 


where n z . is the number of pixels in region i, and x z is the mean 
vector for region i. 

If HSWO segmentation is not implemented recursively, it 
will not be subject to problems with processing window arti- 
facts. 

A description of the HSWO algorithm follows and is 
shown in FIG. 7: 

HSWO Basic Algorithm Description 10: 

1 . Initialize the segmentation by assigning each image pixel a 
region label 11. If a pre- segmentation is provided, label 
each image pixel according to the pre-segmentation. Oth- 
erwise, label each image pixel as a separate region. 

2. Calculate the dissimilarity criterion value between all pairs 
of spatially adjacent regions 12, find the pair of spatially 
adjacent regions with the smallest dissimilarity criterion 
value, and merge that pair of regions 13. 

3. Stop 15 if no more merges are required. Otherwise 14, 
return to step 2. 

Beaulieu and Goldberg did not specify a tie-breaking 
method for step 2, i.e., no procedure is specified for handling 
the case when more than one pair of regions was found to have 
the smallest dissimilarity criterion value. A suitable tie-break- 
ing criterion is as follows: If more than one pair of regions has 
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the smallest dissimilarity criterion, the merge involving the 
region with the largest region label is arbitrarily performed 
first. If this region with the largest region label also has the 
smallest dissimilarity criterion relative to more than one 
region, all of these regions are merged together (effectively, in 5 
parallel). This tie-breaking criteria is arbitrary, and is based 
on convenience of implementation. An additional arbitrary 
convention is that the region with the larger region label is 
merged into the region with the smaller region label (i.e., the 
smaller region label survives the merge). 10 

Beaulieu and Goldbeig show that the HSWO algorithm 
produces the globally optimal segmentation result if the sta- 
tistics at all iterations are independent. Even though the sta- 
tistics at all iterations will generally not be independent for 
natural images, the HSWO approach is still shown to produce 
excellent results. Beaulieu and Goldberg also point out that 
the sequence of partitions generated by this iterative approach 
reflect the hierarchical structure of the imagery data: the 
partitions obtained in the early iterations preserve the small ^ 
details and objects in the image, while the partitions obtained 
in the latter iterations preserve only the most important com- 
ponents of the image. They further note that these hierarchical 
partitions may carry information that may help in identifying 
the objects in the imagery data. ^ 

Beaulieu and Goldberg developed a dissimilarity criterion 
(which they call a stepwise criterion) based on the piecewise 
approximation an image by constant value regions, with the 
value given by the mean value of the region. For the gray scale 
(single spectral band) images they considered, their stepwise 30 
merging criterion, C Ij7 , is: 


Cij 


rii n j 

(fii + rij ) 





While HSWO has advantages over the classical image 
segmentation approaches in terms of optimality, it runs into 
problems when applied to even moderately large data sets. 
This is because HSWO does not allow the merging of non- 40 
spatially connected regions. The large number of regions 
required by HSWO for a complete representation of the seg- 
mentation detail leads to processing inefficiencies in two 
different ways: 

45 

(i) If HSWO is implemented recursively (like basic RHSEG), 
this large number of regions limits min_nregions to be no 
less than chk_nregions/4, where chk_nregions is the num- 
ber of regions required for complete representation, and 
min_nregions is the as defined in the basic RHSEG algo- 5Q 
rithm. 

(ii) If a hierarchical segmentation result is desired, the com- 
putations must be computed for a much larger number of 
region growing iterations. 

HSWO and HSEG without spectral clustering are imprac- 55 
tical for even moderately sized image because of the number 
of regions required to represent a useful level of segmentation 
detail. The optional spectral clustering step in the HSEG 
algorithm is what makes practical use of the HSEG algorithm 
possible for even very large images, with the approximation 60 
of HSEG in the RHSEG recursive implementation. (A recur- 
sive implementation of HSEG is not useful without the spec- 
tral clustering step, due to the high value of min_nregions 
required for even moderately sized data sets.) This is similar 
to the multiple spectral band “square root of the band sum 65 
mean squared error” criterion used in HSEG and explained 
under the parameter dissim_crit below. 


HSEG is an elaboration of HSWO in which an option for 
merging spatially non-adjacent regions is added along with a 
method for selecting a subset of segmentation results to form 
the segmentation hierarchy output 20 shown in FIGS. 8 A and 
8B. The variable spclust_wght is introduced to define the 
relative importance of spectral clustering versus region grow- 
ing. When spclust_wght=0, only merges between spatially 
adjacent regions are permitted. When spclust_wght=l, other 
adjacent and non-adjacent regions are considered equally. 
Other italicized parameters mentioned are defined at the end 
of the specification. 

HSEG Algorithm Description: 

1. Give each image pixel a numeric region label 21. If a 
pre- segmentation is provided, label each image pixel 
according to the pre-segmentation. Otherwise, label each 
image pixel as a separate region. 

2. Calculate the dissimilarity value, dissim_val, between all 
pairs of spatially adjacent regions 22. If spclust_wght>0.0, 
also calculate dissim_val between all pairs of spatially 
non-adjacent regions. 

3. Find the smallest dissim_val between spatially adjacent 

pairs of regions 23 and set nghbr_thresh equal to it. Set 
prev_max_threshold= : 0.0 and 

max_threshold=nghbr_thresh. 

4. If spclust_wght=0.0 24, go to step 7. Otherwise, if 
max_threshold=0.0 25, merge all pairs of regions (spatially 
adjacent or spatially non-adjacent) with dissim_val=0.0 
26. If max_threshold>0.0, merge all pairs of regions (spa- 
tially adjacent or spatially non-adjacent) with 
dis sim_val<spclust_wght * max_thresho Id in order from 
smallest dissim_val to largest dissim_val 27. Update the 
dissim_val’s for spatially adjacent and non-adjacent 
regions between each merge as necessary. (Ties are 
resolved in the same manner used for HSWO.) 

5. If the number of regions remaining is less than or equal to 
the preset value chk_nregions 28, go to step 1 1 . Otherwise, 
update the dissim_vaFs between spatially adjacent pairs of 
regions as necessary 29. 

6. Find the smallest dissim_val between spatially adjacent 

pairs of regions and set nghbrjhresh equal to it 30. If 
nghbr_thresh>max_threshold, set 

prev_max_threshold=max_threshold and then set 
max_threshold=nghbr_thresh. 

7. Merge all pairs of spatially adjacent regions with 
dis sim_val ^max_thresho Id in order from smallest 
dissim_ val to largest dissim_val 31. Update the dissim_ 
val’s for spatially adjacent regions between each merge as 
necessary. (Ties are resolved in the same manner used for 
HSWO.) 

8. If the number of regions remaining is less than or equal to 
the preset value chk_nregions, go to step 11, or if 
spclust_wght=0.0 32, go to step 11. Otherwise, update the 
dissim_val’s between all pairs of spatially adjacent and 
non-adjacent regions 33. 

9. Merge all pairs of spatially adjacent or non-adjacent 
regions with dissim_val^ispclust_wght*max_threshold in 
order from smallest dissim_val to largest dissim_val. (For 
the most part, only spatially non-adjacent merges will 
occur.) 34. Update the dissim_val’s for spatially adjacent 
and non-adjacent regions between each merge as neces- 
sary. (Ties are resolved in the same manner used for 
HSWO.) 

1 0. If the number of regions remaining is less than or equal to 
chk_nregions 35, go to step 11. Otherwise, go to step 6. 

1 1 . If the number of regions remaining is less than or equal to 
conv_nregions 36, save the current region label map to disk 
along with associated region information and STOR If this 



US 7,697,759 B2 


9 


10 


is the first time this step is executed 37, save the current 
region label map to disk along with associated region infor- 
mation 38, and go to step 6. Otherwise, calculate 
tratio=max_threshold/prev_max_threshold 39. (Since this 
instruction cannot be reached first time this step is 
executed, prev_max_threshold is guaranteed to be >0.0.) If 
tratio is greater than the preset threshold convfact 40, save 
the region label map from the previous iteration to disk 
along with associated region information, and go to step 6. 
Otherwise, just go to step 6. 

Note that spclust_wght, chk_nregions, convfact, andconv_ 
nregions are user specified parameters (however, default val- 
ues chk_nregions=64, convfact=1.01, and conv_nregions=2 
are usually satisfactory). If spclust_wght=0.0, HSEG is 
essentially identical to HSWO with convergence checking 
(step 11) for segmentation hierarchy selection. The associ- 
ated region information mentioned in step 1 1 above is the 
region number of pixels list, and optionally includes the 
boundary region map, the region number of boundary pixels 
list, the region mean vector list, the region standard deviation 
list, and the region maximum merging threshold list (for each 
region, the maximum merging threshold at which any merge 
occurred in the formation of the region). 

A recursive formulation of HSEG, called RHSEG 50, is 
shown in FIG. 9. 

1 . Given an input data set X, specify the number of levels of 
recursion required (rnb_levels) and pad the input data set, 
if necessary, so the width and height (or length) of the data 
set can be evenly divided by 2 mb ~ IeveIs ~ 1 51. (A good value 
for mb_levels results in a data section at level =mb_level s 
consisting of roughly 1000 to 4000 data points.) Set 
level=l . 

2. Call rhseg(level,X) algorithm 52. 

3. Execute the HSEG algorithm (as above) 53 using as a 
pre- segmentation the segmentation output by the call to 
rhseg( ) in step 2. (Continue executing HSEG until the 
number of regions reaches chk_nregions 54 and save the 
segmentation results 55 as specified in step 1 1 of the HSEG 
algorithm.) 

An outline of the “switch pixels” rhseg(level,X) algorithm 
60 is shown in FIGS. 10A and 10B and described as follows: 

If level=rnb_levels 61, go to step 3. Otherwise, divide the data 
set into equal subsections 62. For n-spatial dimension data, 
divide the data set into equal sections, halving each spatial 
dimension. For 1 -spatial dimension data, divide into two 
equal sections; for 2 -spatial dimension data, divide into four 
equal sections; and for 3 -spatial dimension data, divide into 
eight equal sections. 

1 . and call rhseg(level+ 1 ,sub_X) 63 for each subsection of the 
data set (represented as sub_X). 

2. After the calls to rhseg( ) for each data set subsection from 
step 1 complete processing, reassemble the data segmen- 
tation results 64. 

3. Execute the HSEG algorithm as described above 65 (using 
the reassembled segmentation results are as the pre-seg- 
mentation when level<rnb_levels), with the following 
modification: Terminate the algorithm when the number of 
regions reaches the preset value min_nregions (if level=l, 
terminate at the greater of min_nregions or chk_nregions). 

4. If level=rnb_levels 66, exit 67. Otherwise, switch the 
region assignment of certain pixels in the following man- 
ner: 

(i) For each region initialize a candidate_region_label_set 
(of C++ data type set) to zero size. Let seam_threshold_ 
factor, region_threshold_factor, and switch_pixels_fac- 
tor be user set parameters (the default value for each 


parameter of 1 .5 works well for many data sets) . Also, let 
max_threshold be the maximum merging threshold 
encountered in step 3 above. For each region, accumu- 
late entries to the candidate_region_label_set in the fol- 
5 lowing two ways: 

a. For each pixel contained in the pairs of rows and 
columns along the seam between the data subsections 
reassembled in step 2 for which the pixel across the 
seam belongs to a different region 68, do the follow - 

10 ing: Calculate the dissimilarity between the pixel and 

its current region (own_region_dissim), and calculate 
the dissimilarity between the pixel and the region of 
the pixel across the seam (other_region_dissim). If 

own_region_dissim>seam_threshold_factor*other_ 

1 5 region_dissim, add the region label of the region of 

the pixel across seam to the candidate_region_la- 
bel_set of this pixel’s region. 

b. Compare each region to every other region 69. If the 
dissimilarity between a pair of regions is less than 

20 region_threshold_factor*max_threshold, add the 

region label of each region to the other region’ s can- 
didate_region_label_set. 

(ii) For each region with a non-zero sized candidate_re- 
gion_label_set, compute the dissimilarity of each pixel 
contained in that region to its current region (own_re- 
gion_dissim) and to each region in the region’s candi- 
date_region_label_set (other_region_dissim) 70. If a 
pixel is found to have own_region_dissim> 
switch_pixels_factor !i! other_region_dissim, switch the 
region label for that pixel to the region with the mini- 
mum other_region_dissim value 71. (In the case of two 
regions in the list having equal other_region_dissim val- 
ues, the region with the smallest region label is chosen.) 
35 (iii) If spclust_wght>0.0 72, exit 73. Otherwise perform 
connected component labeling on the current segmen- 
tation results and relabel accordingly 74. Then execute 
the HSEG algorithm 75 on the data X with the following 
modification: Terminate the algorithm when the number 
40 of regions reaches the preset value min_nregions (if 
level=l, terminate at the greater of min_nregions or 
chk_nregions). Exit 76. 

In step 4 data points on the seam between the data subsec- 
tions are analyzed to find pairs of regions that may contain 
45 pixels that are more similar to the other region. All of the 
regions are also compared to each other directly, and regions 
that are relatively similar to each other are assumed to contain 
pixels that may be more similar to the other region. The first 
approach identifies pairs of regions that occur as neighbors 
50 across the processing window seam that contain pixels that 
are more similar to the other region, no matter how similar the 
region pair is to each other. This approach is most effective at 
eliminating obvious blocky artifacts at the processing win- 
dow seam. The second approach identifies pairs or regions 
55 that are relatively similar to each other, and thus may contain 
pixels that are more similar to the other region, even though 
the pair of regions may not occur as neighbors across a pro- 
cessing window seam. This approach is most effective at 
eliminating more diffuse processing window artifacts. 

60 This is very efficient because it focuses specifically on the 
regions that are involved in the processing window artifacts. 
In addition, performing the pixel switching after performing 
the HSEG algorithm as in step 3 enhances processing effi- 
ciency. Pixel switching (as in step 4) could be performed 
65 before step 3, but then a number of pixels would be switched 
between pairs of regions that would end up getting merged 
together in step 3. 
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FIG. 3 shows the segmentation of the image of FIG. 1 into 
8 regions produced by the “switch pixels” approach for elimi- 
nating processing window artifacts, which is summarized 
above. The program was run with spclust_wght=0. 1 and with 
the other parameters set to program defaults. The processing 5 
window artifacts are clearly eliminated. However, in addition 
to the elimination of the processing window artifacts, the 
nature of the segmentation result shown in FIG. 3 is very 
different from that in FIG. 2. 

There are two main reasons for the difference in the nature to 
of the segmentation results in these two cases. The first reason 
is evidenced by the fact that 96 regions are required in FIG. 2 
to represent a segmentation with equivalent mean-square 
error based dissimilarity to the 8 region segmentation in FIG. 

3. The dissimilarity between the region mean image and the 15 
original image data is 0.444 for FIG. 2 and 0.423 for FIG. 3. 

In the process that produced FIG. 2, merges performed 
within a small processing window cannot be undone as the 
recursion returns to levels with larger processing windows. 
This discourages the merging of regions that originated in 20 
widely separated local processing windows, and results in 
requiring many more regions to describe the image with the 
same fidelity. Because of the efficiency of describing the 
image with the same fidelity in many fewer regions, the 
segmentation result in FIG. 3 is better than the result of FIG. 25 
2 . 

However, the second reason for the difference in the nature 
of the segmentation results is not as positive. The hierarchical 
data segmentation algorithm, HSEG, utilizes the parameter 
spclust_wght to control the relative importance of merges 30 
between spatially adjacent regions and merges between spa- 
tially non-adjacent regions. When spclust_wght=0.0, only 
merges between spatially adjacent regions are allowed. When 
spclust_wght=1.0, merges between spatially adjacent and 
spatially non-adjacent regions are given equally priority. For 35 
values of spclust_wght between 0.0 and 1.0, spatially adja- 
cent merges are given priority over spatially non-adjacent 
merges by a factor of 1 .0/spclust_wght. 

A problem with the “switch pixels” approach to eliminat- 
ing processing artifacts is that it does not have any mechanism 40 
for giving priority to spatial adjacency in its process for 
switching pixels from one region to another. Because of this, 
the segmentation results produced by this approach for all 
values of spclust_wght greater than 0.0 are very similar to the 
results for spclust_wght=1.0. As an example of this, FIG. 4 45 
shows the segmentation of this image into 8 regions produced 
by the “switch pixels” approach described with 
spclust_wght=1.0 and other program parameters set to pro- 
gram defaults. FIG. 3 is much more similar to FIG. 4 than it is 
to FIG. 2. 50 

The current disclosure describes an innovative “split-re- 
merge” approach to eliminating processing artifacts that 
splits certain pixels out from regions they are assigned to and 
remerges them into perhaps another, more similar, region, 
while giving the appropriate priority to spatial adjacency. 55 
FIG. 5 shows the 10 region segmentation the result from 
running RHSEG with the “split-remerge” approach to elimi- 
nating processing window artifacts. The program was run 
with spclust_wght=0. 1 and other program parameters set to 
program defaults. 60 

Comparing the segmentation result from FIG. 5 to the 
results from FIGS. 3 and 4, more coherent regions are obvious 
in FIG. 5, with many fewer small sub-regions. This is the 
effect of spclust_wght=0. 1 in the process that produced FIG. 

5. In this sense, FIG. 5 is more similar to FIG. 2, which also 65 
has more coherent regions with many fewer small sub-re- 
gions. 
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Besides the elimination of processing window artifacts, the 
difference between the segmentation results shown in FIG. 2 
and FIG. 5 is due mainly to a subtle effect of the processing 
window elimination process. As noted earlier, in the process 
that produced FIG. 2, merges performed within a small pro- 
cessing window cannot be undone as the recursion returns to 
levels with larger processing windows. This discourages the 
merging of regions that originated in widely separated local 
processing windows, and results in requiring many more 
regions to describe the image with the same fidelity. However, 
as part of the processing window artifact elimination process, 
certain image pixels are split out from the region they were 
originally assigned to and allowed to become members of a 
different region. This is true for both the “switch pixels” and 
“split-remerge” approaches. This allows for the adjustment of 
the regions to better describe the image in a global sense, 
which makes possible the characterization of the image with 
fewer regions at the same fidelity of representation. The 
essential difference between the “switch pixels” approach 
and the “split-remerge” approach described herein is in how 
the pixels are given a new region assignment. With the “split- 
remerge” approach, appropriate priority is given to spatially 
adjacent region merges in the pixel reassignment, while in the 
“switch pixels” approach spatially adjacent and spatially non- 
adjacent merges are given equal priority in the pixel reassign- 
ment. 

The software innovation described below provides hierar- 
chical segmentations of data that are free of processing win- 
dow artifacts and with maximum flexibility in controlling the 
nature of the segmentation results with the spclust_wght 
parameter. The innovation allows the spclust_wght parameter 
to have full effect over its permissible value range while 
producing results free of processing window artifacts. The 
RHSEG algorithm shown in FIG. 9 is employed with a split- 
remerge version of the rhseg(level,X) algorithm. 

An outline of the split-remerge rhseg(level,X) algorithm 
160 is shown in FIGS. 11A and 11B and described as follows: 

1 . If level=mb_levels 161, go to step 3. Otherwise, divide the 
data set into equal subsections 162 and call rhseg(level+l, 
sub_X) for each subsection of the data set (represented as 
sub_X) 163. 

2. After the calls to rhseg( ) for each data set subsection from 
step 1 complete processing, reassemble the data segmen- 
tation results 164. 

3. Execute the HSEG algorithm (as described in connection 
with FIGS. 8 A and 8B, using the reassembled segmenta- 
tion results are as the pre- segmentation when 
level<mb_levels), with the following modification: Termi- 
nate the algorithm when the number of regions reaches the 
preset value min_nregions (if level=l, terminate at the 
greater of min_nregions or chk_nregions) 165. 

4. If level=mb_levels 166, exit 167. Otherwise, modify the 
region assignment of certain pixels in the following man- 
ner: 

(i) For each region initialize a candidate_region_label_set 
(of C++ data type set) to zero size. Let seam_threshold_ 
factor, region_threshold_factor, and split_pixels_factor 
be user set parameters (the default value for each param- 
eter of 1.5 works well for many data sets). Also, let 
max_threshold be the maximum merging threshold 
encountered in step 3 above. For each region, accumu- 
late entries to the candidate_region_label_set in the fol- 
lowing two ways: 

a. For each pixel contained in the pairs of rows and 
columns along the seam between the data subsections 
reassembled in step 2 for which the pixel across the 
seam belongs to a different region, do the following: 
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Calculate the dissimilarity between the pixel and its 
current region (own_region_dissim), and calculate 
the dissimilarity between the pixel and the region of 
the pixel across the seam (other_region_dissim). If 
own_region_dis sim>seam_threshold_factor *other_ 5 
region_dissim, add the region label of the region of 
the pixel across seam to the candidate_region_la- 
bel_set of this pixel’s region 168. 
b. Compare each region to every other region. If the 
dissimilarity between a pair of regions is less than 10 
region_threshold_factor*max_threshold, add the 
region label of each region to the other region’s can- 
didate_region_label_set 169. 

(ii) For each region with a non-zero sized candidate_re- 
gion_label_set, compute the dissimilarity of each pixel 15 
contained in that region to its current region (own_re- 
gion_dissim) and to each region in the region’s candi- 
date_region_label_set (other_region_dissim) 170. If a 
pixel is found to have own_region_dissim>split_pixels_ 
factor*other_region_dissim, split the pixel out from its 20 
current region and flag as a “split pixel.” 

(iii) If spclust_wght=l .0, for each pixel flagged as a “split 
pixel,” merge it into the region for which it has the 
minimum other_region_dissim value of the regions in 
the candidate_region_label_set of the region it previ- 25 
ously belonged to 173a. (In the case of two regions in the 
list having equal other_region_dissim values, the region 
with the smallest region label is chosen.) Otherwise, 
execute remerge(level,max_threshold,X) on the current 
region labeling (the function remerge( ) is outlined 30 
below). 

(iv) If spclust_wght>0.0 174, exit 175. Otherwise perform 
connected component labeling on the current segmen- 
tation results and relabel accordingly 176. Then execute 
the HSEG algorithm on the data X 177 with the follow- 35 
ing modification: Terminate the algorithm when the 
number of regions reaches the preset value min_nre- 
gions. Exit 178. 

Notice the differences between this split-remerge version 
of RHSEG and the switch pixels version in FIG. 8. In step 40 
4(ii) of the switch pixels version, pixels that have 
own_region_dissim>switch_pixels_factor*other_region_ 
dissim are directly switched from their current region assign- 
ment to the most similar region in the current region’s candi- 
date_region_label_set. In step 4(ii) of the split-remerge ver- 45 
sion, pixels that have own_region_dissim>split_ 
pixels_factor*other_region_dissim are just split out from 
their current region, and flagged as a “split out” pixel (note 
that switch_pixels_factor has been renamed to split_pixels_ 
factor). 50 

In this split-remerge version, a new step 4(iii) is added 
between the step 4(ii) and step 4(iii) of the switch pixel 
version. In the case where spclust_wght=l .0, the pixels that 
were flagged as being “split out” in step 4(ii) are simply 
merged into the region for which it has the minimum other_ 55 
region_dissim value of the regions in the candidate_region_ 
label_set of the region it previously belonged to. This makes 
the combination of steps 4(ii) and 4(iii) for the current version 
equivalent to step 4(ii) of the previous version for 
spclust_wght=1.0. However, for spclust_wght<l .0, the cur- 60 
rent version of RHSEG is very different from the previous 
version. 

For spclust_wght<1.0, the current version executes 
remerge(level,max_threshold,X) on the current region label- 
ing. The remerge( ) function is one of the key innovations, 65 
which directly and through its use of a restructed HSEG 
algorithm introduces contagious regions. This recursive 


remerge(level,max_threshold,X) function 80 is outlined 

below and illustrated in FIG. 12: 

1. If level<rnb_levels 81, divide the data set into equal sub- 
sections (as in rhseg( ) and call remerge(level+l, 
max_threshold,sub_X) 83 for each subsection of the data 
set (represented as sub_X) 82, and go to step 2. Otherwise, 
let onregions=current number of regions, which are labeled 
1 to onregions 84. Give each split out pixel a numeric 
region label, starting at onregions+1, and flag them as a 
new_region 85. Give each new region a candidate_region_ 
label_set equal to the candidate_region_label_set of the 
region it was split out from and add to that set the region 
label of the region it was split out from 86. 

2. After the calls to remerge( ) for each data set subsection 
from step 1 complete processing, reassemble the data seg- 
mentation results 87 (for all but the first subsection, the 
regions flagged as new_region must be renumbered to not 
duplicate labels of the new regions from the previous sub- 
sections). 

3. If level>the value of level at the initial call to remerge 88, 
flag as contagious each region flagged as a new_region that 
touches an interior processing window. Then, execute a 
restricted version of the HSEG algorithm as described 
below 89ci (if level<rnb_levels, use the reassembled seg- 
mentation results a pre-segmentation), terminating when 
the number of regions equals converge_nregions= 
onregions+min_nregions 89 b. Exit. 

An outline of restricted HSEG 90 follows and is illustrated 

in FIGS. 13A and 13B: 

1. For all regions flagged as a new_region, calculate the 
dissimilarity value, dissim_val, between it and all spatially 
adjacent regions (including regions not flagged as a 
new_region) 91. 

2. Out of the regions flagged as a new_region, find the region 
with smallest dissim_val 92 (ties are resolved in the same 
manner used for HSWO). If dissim_val>max_threshold 
93, go to step 5. Otherwise, call this region best_region and 
its most similar neighboring region merge_region 94. 

3. If best_region is flagged as contagious 95, drop best_region 
from further consideration of spatially adjacent merging 
(for this call to restricted HSEG) 96, and, if merge_region 
is flagged as a new_region, flag it as contagious. Go to step 
2 . 

4. Merge best_region into merge_region 97 (if merge_region 
is flagged as a new_region, it retains its new_region flag, 
and the candidate_region_label_set’s of the two regions 
are combined). If the number of 
regions>converge_nregions 98, update as necessary the 
dissim_val’s of regions flagged as a new_region and their 
spatially adjacent regions 99, and go to step 2. Otherwise, 
exit. 

5. Returning the regions previously dropped from merging 
back into consideration, calculate the dissimilarity value, 
dissim_val, between the regions flagged as a new_region 
and all spatially adjacent regions (including regions not 
flagged as a new_region). Also calculate the dissimilarity 
value, dissim_val, between the regions flagged as a 
new_region and all regions in their candidate_region_la- 
bel_set 100. Let last_ngbhr_thresh be the threshold of the 
last merge that occurred in step 4 above. 

6. Out of the regions flagged as a new_region, find the region 

with smallest dissim_val 101 (ties are resolved in the same 
manner used for HSWO). If 

dissim_val*spclust_wght>last_nghbr_threshold 102, exit. 
Otherwise, call this region best_region and its most similar 
region merge_region 103. 
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7. Merge best_region into merge_region 104 (if merge_re- 
gion is flagged as a new_region, it retains its new_region 
flag, and the candidate_region_label_sets’s of the two 
regions are combined). If the number of 
regions<=converge_nregions 105, exit. Otherwise, update 
as necessary the dissim_val’s of regions flagged as a 
new_region and their spatially adjacent regions and all 
regions in their candidate_region_label_set 106. Go to step 
6 . 

The remerge( ) function is called recursively from the 
recursive level it was initially called from in step 4(iii) of the 
rhseg( ) function. At recursive level mb_levels specially 
flagged single pixel regions are formed from the pixels that 
were split out from their region assignment in the current 
processing window. These new regions have their new_re- 
gion flag set true, and any of these new regions that are located 
on an interior processing window boundary have their conta- 
gious flag set true. In addition, each of these new regions have 
their candidate_region_label_set set equal to the candidate_ 
region_label_set of the region it previously belonged to aug- 
mented by the label of the region it previously belonged to (to 
allow for the possibility of merging back into the previous 
region). 

In the restricted version of HSEG, nearest neighbor merges 
are considered first up to the maximum merge threshold pre- 
viously encountered in the last call to the unrestricted version 
of HSEG. These nearest neighbor merges are constrained by 
delaying merges involving regions touching the interior pro- 
cessing window boundaries until those regions no longer are 
on a processing window boundary or the recursive level is 
reached from which the remerge( ) function was initially 
called from. 

If the number of regions converge_nregions is not reached 
during the processing of the nearest neighbor merges, more 
general merges are considered. These merges are not con- 
strained by the contagious flag, and consider possible merges 
in each region’ s candidate_region_label_set in addition to the 
neighboring regions. In this case the merges are constrained 
by last_nghbr_threshold * spclust_wght as well as converge_ 
nregions. 

The contagious region aspect of the restricted version of 
HSEG is a very important part of this new innovation in the 
processing window artifact elimination process. Without it, 
processing window artifacts are reintroduced, as demon- 
strated by the results displayed in FIG. 6. 

When the contagious region concept is incorrectly imple- 
mented, it will be less computationally efficient than the 
current approach. This is because when the contagious region 
idea is applied to the region growing process starting at single 
pixel regions over the entire processing window, considering 
contagious regions at the interior processing window bound- 
aries invariably causes merging to stop at a number of regions 
much higher the usual value of min_nregions utilized in the 
current version of RHSEG. Larger values of min_nregions 
have a deleterious effect on processing time. However, the 
contagious region approach does not lead to these problems in 
the current implementation since it is applied only to a subset 
of regions corresponding to a very limited subset of image 
pixels, and since after nearest neighbor merges are stopped by 
all new regions becoming contagious, additional general 
merges may be performed. These additional merges do not 
lead to artifacts because for the most part the merges involve 
merges of the new regions into globally defined “old” regions, 
and because a large number of nearest neighbor merges have 
already occurred. 

As shown in Table 1, the processing times required by 
HSEG and RHSEG for spclust_wght>0.0 are generally 
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longer than required with spclust_wght=0.0. The contrast in 
processing times is most pronounced for HSEG (rnb_lev- 
els=l) as the image size grows. Even for a very moderate 
image size, 256x256 pixels, the processing time required by 
5 HSEG for spclust_wght>0.0 is about 235 times longer than 
that required by HSEG for spclust_wght=0.0. This is because 
the number of regions that need to be compared in the early 
stages of HSEG is on the order of N 2 when spclust_wght>0.0, 
while the number of regions that need to be compared in the 
to early stages of HSEG is on the order of N*n when 
spclust_wght=0.0, where N is the number of pixels in the 
image and n is the size of a pixel neighborhood (eight neigh- 
bors, by default). This increase in processing time is much 
less pronounced for RHSEG using the default values for 
15 mb_levels. In this case, for the image sizes tested, the pro- 
cessing time for the parallel version of RHSEG with 
spclust_wght>0.0 is up to about 4 to 5 times longer than that 
for RHSEG with spclust_wght=0.0. This is because the num- 
ber of regions that need to be compared in the early stages of 
20 RHSEG is limited, by default, to no more than 2048 regions. 

Functional Operation of the Algorithm in Software and 
Parameters 

Provided in this section is a user’ s guide-like description of 
25 the parameters used in the HSEG and RHSEG algorithms 
(HSEG is simply RHSEG with rnb_levels=l) embedded in a 
software program. This implementation assumes the input 
data is one or two -spatial dimension single-band, multispec- 
tral of hyperspectral image (or image-like) data. Selection 
30 between the single processor version and the multiprocessor 
(parallel) version is made through use of a compiler prepro- 
cessing directive. 

The single processor version of RHSEG is run from the 
command line with the command: 

35 rhseg parameter_file_name 

where parameter_file_name is the name of the input param- 
eter file (for contents, see below). The parallel (multiproces- 
sor) version of RHSEG is run from the command line with the 
command: 

40 rhseg parameter_file_name inb_levels onb_levels 

where parameter_file_name is the name of the input param- 
eter file (for contents, see below). The other command line 
parameters for the parallel version are: 


inb_levels (short unsigned int) Data input recursive level. 

(1 = inb_levels = rnb_levels) 


50 

(See below for the description of mb_levels.) This is the 
recursive level at which the input data is input. Process- 
ing at recursive level s>inb_levels is performed sequen- 
tially. This parameter is used to optimize the parallel 
55 processing efficiency. The optimal value depends on the 
parallel processing hardware and the size of the data 
being processed. 

The value of inb_levels determines the number of proces- 
sors (numprocs) utilized in the parallel version. For two- 
60 spatial dimension data as input, numprocs=4 m£,_/eveZy-1 
For 1 -dimensional data, numproc s=2 zw ^ e v ’ eZs “ 1 , and for 
3-dimensional data, numprocs=8^ eve ^- 1 . 

At program initialization, the input data is parceled out to 
the numprocs processors, and each processor indepen- 
65 dently processes each of the numprocs data sections. 
After each process finishes with its section of data, it 
transfers the results to the appropriate processor work- 
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ing at recursive level inb_levels-l, etc. until the final 
stage of processing at recursive level 1 . 


onb_levels (short unsigned int) Data output recursive level. 

(1 onb_levels = inb_levels) 


5 


The lowest recursive level at which the input data and other 
pixel oriented data (such as the region label map) is 10 
maintained. As such, it is the recursive level from which 
the region label map and other pixel-oriented outputs are 
output from. This parameter is used to optimize the 
parallel processing efficiency. The optimal value 
depends on the parallel processing hardware and the size 
of the data being processed. 
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Pat. No. 6,895,1 15 to reuse processes at different levels 
of recursion. In the illustrated case, onb_levels is 2 
(2x2=4 data sections on 4 tasks), inb_levels is 3 (4x4=16 
data sections on 16 tasks), mb_levels is 5 (16x16=256 
data sections on 16 tasks), and the input image has 512 
columns and 512 rows. At program initialization, the 
input data is split into 16 sections and distributed to the 
parallel tasks at recursive level inb_levels. From that 
point, each task recursively divides its data in order to 
sequentially process its data from recursive levels 
mb_levels back up to inb_levels. Once the recursive 
processing is completed up to recursive level inb_levels, 
the input and region label map results are sent to tasks 
running at recursive level inb_levels-l , etc., until recur- 
sive level onb_levels is reached. Processing then contin- 


TABLE 1 


Processing times for HSEG and RHSEG for variations in number of processors, 

the number of levels of recursion (mb levels, where rnb levels = 1 means 

no recursion), and spclust wght. Blank entries signify that the processing 

would take more than 1 hour. 

Timings were performed with 2.4 GHz processors where 1 Gbyte of memory per CPU 
was available for the single processor runs, and 0.5 Gbyte of memory per CPU 
was available for the multiple processor runs. 

Processing Times (minutes :seconds) 


Number of spclust wght = spclust wght = spclust wght = 

Image Size Processors rnb levels 0.0 0.1 1.0 


0064 x 0064 

1 

1 

0:01 

0:04 

0:05 

0128 x 0128 

1 

1 

0:02 

1:14 

1:14 

0256 x 0256 

1 

1 

0:05 

19:28 

19:36 

0512 x 0512 

1 

1 

0:20 

— 

— 

1024 x 1024 

1 

1 

1:24 

— 

— 

0064 x 0064 

1 

2* 

0:01 

0:04 

0:03 

0128 x 0128 

1 

3* 

0:02 

0:18 

0:13 

0256 x 0256 

1 

4* 

0:08 

1:47 

0:56 

0512 x 0512 

1 

5* 

0:51 

8:51 

3:54 

1024 x 1024 

1 

6* 

5:29 

40:11 

16:03 

2048 x 2048 

1 

7* 

43:57 

— 

— 

0064 x 0064 

4 

2 * 

0:01 

0:02 

0:02 

0128 x 0128 

16 

3* 

<0:01 

0:04 

0:03 

0256 x 0256 

64 

4* 

0:02 

0:06 

0:05 

0512x0512 

256 

5* 

0:02 

0:11 

0:09 

1024 x 1024 

256 

6* 

0:05 

0:23 

0:19 

2048 x 2048 

256 

7* 

0:16 

1:18 

0:54 

4096 x 4096 

256 

8* 

1:08 

4:35 

2:45 

6912 x 6528 

256 

9* 

3:36 

8:23 

4:15 


* Default value for mb_levels for this image size. 


When RHSEG reaches processing at recursive level 
onb_levels, the input data, region label map, and any 
other pixel oriented information is maintained in the 
processors operating at recursive level onb_levels, and 
not transmitted up the recursive processing chain. When 
information relating to the input data or region label map 
is required at processing levels<onb_levels, this infor- 
mation is passed from the processors active at the 
onb_levels recursive level through interprocessor com- 
munication up to the appropriate process at the current 
recursive level. 

Setting onb_level>l distributes the pixel -oriented data 
over 4 ? nb - 2evels - 1 processors, and enables the processing 
of very large images by reducing the computer memory 
demands per processor. 

See FIG. 14 for a graphical depiction of the relationship 
between recursive levels rnbjevels, inb_levels and 
onb_levels. The process numbering for this parallel pro- 
cessing scheme is modified from that introduced in U.S. 


ues through recursive level 1, but the input data and 
50 region label map results remain at recursive level 
onb_levels. Note that process 0 is active at all levels of 
recursion, and processes divisible by 4 (corresponding 
to the onb_levels level of recursion) are active at all 
levels of recursion, but is otherwise unchanged. 

NOTE: The parameters inb_levels and onb_levels are used 
only in the parallel version and are not defined or used in the 
single processor version. 

The following required parameters specify and describe 
60 the input data: 


indata (string) Input image data file name 

65 

The input image data file from which a hierarchical image 
segmentation is to be produced. This image data file is 
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assumed to be a headerless binary two-spatial dimension 
image or image-like data file in band sequential format. 
The number of columns, rows, spectral bands and the 
data type are specified by other required parameters (see 
below). Data types “unsigned char (byte)” and “short 5 
unsigned inf ’ are supported. 


ncols (int) Number of columns in input image data 10 

(0 < ncols < 65535) 

mows (int) Number of rows in input image data 

(0 < mows < 65535) 

nbands (short unsigned Number of spectral bands in input image data 
int) 

(0 < nbands <65535) ^ 

dtype (short unsigned Data type of input image data 
int) 

dtype = 8 designates “unsigned char (byte)” 
dtype =16 designates “short unsigned int” 

(otherwise undefined) 


20 

The following required parameters specify output files: 


rlblmap (string) 


Output region label map data file name 25 


The region label map at the finest level of segmentation 
detail (hierarchical level 0). Together with regmerges 
(see below), this forms the main output of RHSEG. 
Region label values of “0” correspond to invalid input 30 
data values in the input image data. Valid region label 
values range from 1 through 65535. The data is of data 
type “short unsigned int.” 


35 
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The output parameter file contains (in ASCII form) all the 
output parameters from RHSEG. This parameter file is 
formatted in the same way as the input parameter file for 
RHSEG and contains most of the same parameters. 
Additional parameters are the number of hierarchical 
segmentation levels (nb_levels) in the hierarchical seg- 
mentation output and the number of regions 
(level0_nregions) in the hierarchical segmentation with 
the finest segmentation detail. These additional param- 
eter values are required to interpret the rnpixlist, reg- 
merges, rmeanlist, rthreshlist and boundary_npix output 
files (see below). 

The following parameters specify recommended, but 
optional, output files (no defaults): 


rmeanlist (string) Output region mean list file name 

(default = {none}) 


The region mean list file is an optional output of RHSEG. 
This list consists of the region mean value (of data type 
“double”) of each region stored as rows of values and 
groups of rows, with the column location (with counting 
starting at 1) corresponding to the region label value, the 
row location (in each row group) corresponding the 
spectral band and row group corresponding to the seg- 
mentation hierarchy level (with counting starting at 0). 


rthreshlist 


(string) 


Output region maximum merge thr eshold 
list file name 

(default = {none}) 


regmerges (string) Output region label map merges list file name 


The region merges list file consists of the renumberings of 
the region label map required to obtain the region label 40 
map for the second most detailed level (hierarchical 
level 1 ) through the coarsest (last) level of the segmen- 
tation hierarchy from rlblmap (see above). The data type 
is “short unsigned int.” The data is stored as rows of 
values, with the column location (with counting starting 45 
at 1 ) corresponding to the region label value in rlblmap 
(the region label map at the finest level of detail of the 
segmentation hierarchy) and the row location corre- 
sponding to the segmentation hierarchy level (the \ th row 
is the renumberings required to obtain the (1+ 1 ) th level of 50 
the segmentation hierarchy). 


rnpixlist (string) Output region number of pixels list file name 


The region number of pixels list consists of the number of 
pixels (of data type “unsigned int”) in each region stored 
as rows of values, with the column location (with count- 
ing starting at 1 ) corresponding to the region label value 60 
and the row location corresponding to the segmentation 
hierarchy level (with counting starting at 0). 


The region maximum merge threshold list file an optional 
output of RHSEG. This list consists of the maximum 
merge threshold encountered in all merges involving 
each region. The values (of data type “double”) are 
stored as rows of values, with the column location (with 
counting starting at 1 ) corresponding to the region label 
value and the row location corresponding to the segmen- 
tation hierarchy level (with counting starting at 0). 


boundary_npix (string) Output region number of boundary pixels list 
file name 

(default = {none}) 


The region number of boundary pixels list is an optional 
output of RHSEG. This list consists of the number of 
boundary pixels in each region (of data type “unsigned 
int”) stored as rows of values, with the column location 
(with counting starting at 1) corresponding to the region 
label value and the row location corresponding to the 
segmentation hierarchy level (with counting starting at 
0 ). 


boundary_map (string) Output hierarchical boundary map file name 

(default = {none}) 


oparam (string) Output parameter file name 65 The hierarchical boundary map is an optional output of 

RHSEG. The data values of this map are (of type 

unsigned char (byte)) 



21 


US 7,697,759 B2 


reg_std_dev (string) Output region standard deviation value list file 
name 

(default = {none}, 
ignored if spatial_wglit = 0.0) 


The region standard deviation value list is an optional 
output of RHSEG. This list consists of the region’s stan- 10 
dard deviation value (of data type “double”) stored as 
rows of values, with the column location (with counting 
starting at 1 ) corresponding to the region label value and 
the row location corresponding to the segmentation hier- 
archy level (with counting starting at 0). 15 

The following parameters specify optional input files and 
associated parameters (with defaults, if any): 


22 

The following are optional parameters are recommended for 
variation by all users (defaults provided): 


spclust_wght (float) Relative importance of spectral clustering versus 
region growing (spatial_wght ^ 0.0, 
default = 1 .0) 

dissim_crit (short unsigned int)Dissimilarity 

criterion“l -Norm,” 

1 . “ 2 -Norm,” 

2 . “Infinity Norm,” 

3 . “(undefined),” 

4 . “(undefined),” 

5 . “Square Root of Band Sum Mean Squared 
Error,” 

6. “Square Root of Band Maximum Mean 
Squared Error,” 

7 . “(undefined),” 

8. “SAR Speckle Noise Criterion.” 

(default: 6 “Square Root of Band Sum Mean 
Squared Error”) 
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mask (string) Input data mask file name (default = {none}) 


The optional input data mask must match the input image 
data in number of columns and rows. Even if the input 25 
image data has more than one spectral band, the input 
data mask need only have one spectral band. If the input 
data mask has more than one spectral band, only the first 
spectral band is used and is assumed to apply to all 
spectral bands for the input image data. If the data value 30 
of the input data mask is not equal to mask_value (see 
the next parameter definition), the corresponding value 
of the input image data object is taken to be a valid data 
value. If the data value of the input data mask object is 
equal to mask_value, the corresponding value of the 35 
input image data obj ect is taken to be invalid and a region 
label of “0” is assigned to that spatial location in the 
output region label map data. The input data mask data 
type is assumed to be “unsigned char.” 


Criterion for evaluating the dissimilarity of one region from 
another. 

Dissimilarity criteria 1 , 2 and 3 are based on vector norms. 
The 1-Norm of the difference between the region mean vec- 
tors, u z and u j7 of regions X z - and X 7 , each with B spectral 
bands, is: 


Ik -k/IIi = ^ \&b -Pjbl 

b= 1 


where \\, ib and \i Jb are the mean values for regions i and j , 
respectively, in spectral band b. The dissimilarity function for 
regions X. and X., based on the vector 1-Norm, is given by: 

d\-Norm(Xi>Xj) = \\ U r U j\\l- ( 4 b) 


The vector 2-Norm of the difference between the region mean 

mask_value If input data mask file is provided, this is the value vectors, U z and U^, of regions X z and is. 

(short unsigned int) in the mask file that designates bad data. Otherwise 
this is the value in the input data that designates 
bad data. (If mask file provided, default = 0 . 45 

Otherwise no default.) 

rlblmap_in Input region label map file name. Data type must be 

(string) short unsigned int (default = {none}) 


^ z (ftib ftjb) 


The optional region label map must match the input image 50 
data in number of columns and rows. If provided, the 
image segmentation is initialized according to the input 
region label map instead of the default of each pixel as a 
separate region. Wherever a region label of “0” is given 
by the input region label map, the region labeling is 55 
assumed to be unknown and the region label map is 
initialized to one-pixel regions at those locations (except 
see rlblmap_mask_value below). 

60 

rlblmap_mask_value If input region label map file is provided, this is 
(short unsigned int) the value in the region label map file that 

designates bad data. NOTE: The output region 

label map file uses the value zero (0) to designate 

bad data, (default = {none}) ^ 


The dissimilarity function for regions X z and X ; , based on the 
vector 2 -Norm, is given by: 

di^QWjHurUjh- ( 5 b) 

The vector 00 -Norm of the difference between the region 

mean vectors, u z . and u^., of regions X z and X^. is: 

lk- M Joo=max( I, b=l, 2 , . . . , B), (6a) 

The dissimilarity function for regions X z and X 7 , based on the 
vector 00 -Norm, is given by: 

d^^x^uru^. (6b) 

Dissimilarity criteria 6 and 7 are based on minimizing the 
increase of mean squared error between the region mean 
image and the original image data (9) (10) (11). The sample 
estimate of the mean squared error for the segmentation of 
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band b of the image X into R disjoint subsets X 1; X 2 , . . . , X^ Since 

is given by: 


l R 

MSE„(X) = — £ MSE b (Xi), 


n-Mb + tjfijb 


MSE b (X i )= y ( Xpb-Kb f 


an alternate form for Equation (9a) is: 


is the mean squared error contribution for band b from see- A v , „ ,, ,2 , „ ,,, ,, ,2 

tt n ± r * - -1 . A MSt b {Xi, Xj) = niiliit - fi Ub ) + nj{fi jb - fi ijb ) = 

ment X z . Here, x is a pixel vector (in this case, a pixel vector 
in data subset X J, and % pb is the image data value for the b th n ^lb ~ ^iEibEjb + m^l + njtf jb - 2njfi jb fi ijb + = 

spectral band of the pixel vector, x . A dissimilarity function 2 2 _ 9r 

based on a measure of the increase in mean squared error due 15 n,Mlb + Uj ^ jb n '^ lb + Hj ^ jb ^ ljb + 

to the merge of regions X, and X, is given by: , )L +B; +B 


d B SMSE(Xi, Xj) = y A MSEtiXi, Xj), 


(«/ + nj)$ jb — — — 7 [(«/ + nj)(riifif b + njfi 2 jb ) - 
yli + n j) 

2(nifi ib +njHj b ) 2 + (iiiftt, +njju jb ) 2 ] = 


- [(«; + nj)(ninf b +njH 2 Jb ) - (mpa, + njfi jb ) 2 ] = 


A MSE b ( X t , X 7 ) = MSE b (X; U Xj) - MSE b {Xi) - MSE b (Xj). (8b) 


7 ( n iEib + n i n jEj b + n i n jlil + 


njltfb ~ n fA ~ 2ninjii ib iij b - n]ti 2 jb ) = 


BSMSE refers to “band sum MSE .” Instead of summing over 
the bands in (7a) one could take the maximum over the 
spectral bands, resulting in a “band maximum MSE:” 

dsMMSEiXi, Xj.)=max{A MSE b {X b Xj), b=l, 2 ...,B}. (8c) 

Using (7b) and exchanging the order of summation, (8b) 
can be manipulated to produce an efficient dissimilarity func- 
tion based on aggregated region features: 


-?w>) = toT? )</* - ma) ■ 


(8c| Combining Equations (8a) and (9c), 


dBSMSE(X» Xj) = A (ftb-Pjb) ■ 

V 1 b= 1 


A MSE b (Xi, Xj) = 


Y \-Xpb - Eijb] 2 - Y bCpb-Eibf- Yj t Xpb-t*jb ] 2 = 


Similarly combining Equations (8c) and (9c), 


^ [(Xpb Eijb) (Xpb Eib) ] + 


dBMMSEiXi, Xj) = -m.ax{(jt ib -ji Jb ) 2 \ b- 1, 2, ... , B). 

J («/ + rij) J 


\kXpb Eijb ) (Xpb fljb) ] — 

XpGXj } 

Y ^X 2 pb - 2 XpbEijb + tfjb - xlb + 2 XpblHb ~ fill + ' 

XpGXi 

Y t*lb - 2 X P bEijb + Eyb ~ X 2 P b + 2 X P bftjb - fi 2 jb \ 


The dimensionality of the d bsmse an d the d BM mse dissimi- 
45 larity criteria is equal to the square of the dimensionality of 
the image pixel values, while the dimensionality of the vector 
norm based dissimilarity criteria is equal to the dimensional- 
ity of the image pixel values. To keep the dissimilarity criteria 
dimensionalities consistent, HSEG uses the square root of 
these dissimilarity criteria. The “Square Root of Band Sum 
Mean Squared Error” criterion is: 


' -2 Eijb Y Xpb + n i(4jb + 2 Hib Y Xpb ~ nu 


~ 2 Eijb Y Xpb+njH 2 j b +2njb Y Xpb-njft 2 b 


4 / S 2 MSE(X i ,Xj)= ^^Y&ib-fljb) 2 , 

W + n j) b=l 


- 2niftibHijb + nmfjb + 2nnif b - nml 'j 
-2njjij b jij b + nji4 jb + 2rijii) b - njit 2 jb J 


tiiinl ~ ^EibEijb + l4jb) + n M)b ~ 2 E jbEijb + l4jb) = 


and the “Square Root of Band Sum Maximum Squared Error” 


H ( Eib ~ Eijb ) +nj(n jb - jt ijb ) , 


where \i iJb is the mean value for the b th spectral band of the 
mean vector, u y , of region represented by X^X^UX,-. 


cntenon is: 


dBMMSEiXi, Xj) = -^) z : b = 1, 2, ... , £}j 
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Dissimilarity criterion 9 is based on the “SAR Speckle 
Noise Criterion” from Beaulieu. The criterion is: 


dsAR{X h Xj) = [ninjifii + «/)] 1/2 ^ 


Iftjb ~ ftjbl 
+njiij b y 


(12) 5 


NOTE: Other dissimilarity criterion can be included as addi- io 
tional options without changing the nature of the RHSEG 
implementation. 


chk__nregions Number of regions at which convergence 

(short unsigned int) factor checking is initiated 

(2 =1 chk_nregions < 65535, default = 64) 
convfact (float) Convergence factor 

(1 convfact = 100.0, default =1.01) 

20 


The following optional parameters may need to be modified 
depending on your operating system (defaults provided): 
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byteswap in Flag specifying byteswapping for input data 

(short unsigned int) (1 (true) or 0 (false), default = 0) 

byteswap out Flag specifying byteswapping for output data 

(short unsigned int) (1 (true) or 0 (false), default = 0) 

30 


The default values should be used for the following optional 
parameters, except in special circumstances (defaults pro- 
vided): 

35 


spatial wght (float) Weight for standard deviation spatial feature 

(spatial_wght = 0.0, default = 0.0) 


Setting spatial_wght=l .0 weights the spatial feature equally 
with the spectral band features, spatial_wght<l .0 weights the 
spatial feature less and spatial_wght>l weights the spatial 
feature more. If D is the dissimilarity function value before 
combination with the spatial feature value, the combined 45 
dissimilarity function value (comparing regions i and j), D c , 
is: 


D c =D+spatial_wght* I sf-sfj I (13) 

50 

where sf . and st) are the spatial feature values for regions i and 
j, respectively. 

The spatial feature employed is the spectral band maxi- 
mum region standard deviation. For regions consisting of 2 or 
more pixels, the region standard deviation for spectral band b 55 
of region i is: 
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where n z is the number of pixels in the region and \.i lh is the 
region mean for spectral band b of region i: 

Hib = ~ X Xpb ' 

1 Xp^Xi 


The spatial feature value for region i is then defined as: 

sfrVi = B*max{o ib :b=\, 2, . . . , B} (15) 

where B is the number of spectral bands. 

The region standard deviation is not defined for regions 
consisting of only one pixel. Further, the region standard 
deviation as calculated by equation (14) can only be consid- 
ered a rough estimate for small regions (say, regions less than 
9 pixels in size). Thus, if one of the regions being compared 
consists of less than 9 pixels, the spatial_wght factor is modi- 
fied by a std_dev_factor as follows: 

spatial_wght'=std_dev_factor*spatial_wght, (1 6a) 


where 


std_dev_factor=(min_npix- 1 . 0)/8 .0, (16b) 

and min_npix is the number of pixels in the smallest of the 
two regions being compared. Note that for min_npix=l, std_ 
dev_factor=0.0. Thus, std_dev_factor serves to gradually 
phase in the standard deviation spatial feature as the regions 
get larger. 


Neighbor connectivity type 

conn_type (short unsigned int) (2-spatial dimensions) 


1. "Four Nearest Neighbors,” 

2. "Eight Nearest Neighbors,” 

3. "Twelve Nearest Neighbors ,” 

4. “Twenty Nearest Neighbors,” 

5. “Twenty-Four Nearest Neighbors,” 
(default: 2 “Eight Nearest Neighbors”) 


based on the following neighborhood chart, where the focal 
pixel is marked “X”: 
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Using this chart, N Nearest Neighbors include pixels 1, 
2, . . . N. 

For 1 -spatial dimension data, cases 1 and 2 degenerate to 
“Two Nearest Neighbors” and cases 3, 4 and 5 degenerate to 
“Four Nearest Neighbors,” according to the following neigh- 
borhood chart: 

II i 1 O i 4 i 


o- ib = 



(X P b -Mb) 2 


1 

ni - 1 


Z (Xpb) 2 -m(Mb) 2 

,GX; 


(14) 


normind (short unsigned int) Image normalization type 

1. "No Normalization,” 

2. "Normalize Across Bands,” 

3. “Normalize Bands Separately” 

(default: 2 “Normalize Across 
Bands”) 
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Let % pb be the original value for the p \ th pixel (out of N pixels) 
in the b th band (out of B bands). The sample mean and sample 
variance of the b th band are 


28 

spclust_wght>0.0. If spclust_wght=0.0, spclust_start=0. 
Otherwi se, spclust_start=2 * min_nregions . 


5 

l N { N ( 17 ) 

«- = jjTiXi* aada i = Ocpb-n,) 2 , 

P = 1 P = 1 

10 

respectively. The following transformation of the data, % pb , 
will produce image data, with mean, M, and variance, 2: 


%pb = C X P b - Pb)\ + M = ^ Xpb + M bi where 


( 18 ) 


15 


V = — and M’ b = M - Y — . 
4-J cn, b <Th 


( 19 ) 


20 


Usually, the data is normalized so that 2 2 (=Z)=1, and M=0. 

As written above, the normalization is applied to each 
spectral band separately. It can also be defined to apply 2 s 
equally across all spectral bands. For this case, use 
a=max{a^:b=l, 2, . . . , B} in (18) and (19). However, this 
type of normalization will produce the same hierarchical 
segmentation result as no normalization at all: the dissim_ 
val’s will change, but the tratio values will end up being the 30 
same (see step 1 1 of the “HSEG Basic Algorithm Descrip- 
tion” described in connection with FIG. 8. 


gdissim_crit (short unsigned int) Global dissimilarity criterion 

1. "1-Norm,” 

2. "2-Norm,” 

3. "Infinity Norm,” 

4. “(undefined),” 

5. “(undefined),” 

6. “Square Root of Band Sum Mean 
Squared Error,” 

7. “Square Root of Band Maximum 
Mean Squared Error,” 

8. “(undefined),” 

9. “SAR Speckle Noise Criterion.” 
(default = {none}) 


Criterion for evaluating the quality of the image segmenta- 
tions based on the global dissimilarity of the region mean 
image versus the original image data. 

The global dissimilarity criteria 1, 2 and 3 are based on 
vector norms. The global dissimilarity function, based on the 
vector 1 -Norm, for the R region segmentation of the N pixel 
data set X is given by: 




(20) 


rnb levels (short unsigned int) Number of recursive levels 

(IS rnb levels < 255, default 

calculated) 


The global dissimilarity function, based on the vector 
2-Norm, for the R region segmentation of the N pixel data set 
X is given by: 


The number of recursive levels. The default is calculated such 
that the number of data points in the subsections of data 
processed at recursion level rnb_levels are no more that 2048 
data points. The number of columns and rows at recursion 
level mbjevels is sub jcols=ncols/2'"" and 

sub_nrows=nrows/2'” 6 - fevefa - 1 . 


Z ii x p~ u ^- 

i = l Xp G X} 


45 The global dissimilarity function, based on the vector 
oo-Norm, for the R region segmentation of the N pixel data set 
X is given by: 


min_nregions (unsigned int) Number of regions for convergence in 
recursive stages 

(0 < min_nregions < 65535, default 
calculated) 


„(*) = 


— v y \\x p -ui\\ 


N id 




( 22 ) 


If not specified, the default is calculated to be 
min_nregions=sub_ncols*sub_nrows/4 (for sub_ncols and 55 
sub_nrow see the rnbjevels parameter). 


spclust_start (short unsigned int) Number of regions at and below 

which spectral clustering is utilized. 
Otherwise, the spectral clustering step 
is skipped. (0 ^ spectral_start < 
65535, default calculated) 


65 

The default for spclust_start is calculated such that spectral 
clustering is utilized only part of the time when 


The global dissimilarity criteria 6 and 7 are based on the 
square root of the mean squared error between the region 
mean image and the original image data. The global dissimi- 
larity criterion “Square Root of Band Sum Mean Squared 
Error” is: 


D BSMSe(X) - 


1 

(ATM) 


Z Z Z OCrb-W,? 


Xp e X; b= 1 


( 23 ) 



US 7,697,759 B2 


29 

The global dissimilarity criterion “Square Root of Band 
Maximum Mean Squared Error” is: 


*(*) = 


(24) 5 


1 R 

— Yj ma M(Xpb-Vib) 2 -- b = 1, 2, ■■■ ,B} 


{N- 1) 


i=l x p GXi 


10 


Dissimilarity criterion 9 is based on the “SAR Speckle 
Noise Criterion” from M. Beaulieu. The global criterion is: 


1 R B 

d sar(X)=~Y j Yj ^\Xpb-Vibl 


/=! x p &Xi b= 1 


15 


(25) 


20 


min_npixels_pct (float) 


Number of pixels for a region as a percentage 
of the total number of pixels in a processing 
window below which merging is accelerated. 
For regions smaller than this minimum number 
of pixels, the dissimilarity function is 
adjusted to favor merging. 

(0.0 < min_npixels_pct < 100.0, 
default = 0.0) 


30 


This parameter is used to calculate min_pixels as follows: 

m in jipixel s= [ np ixe Is * min_np ixels _p ct/\ 0 0 . 0 J , (26a) 

where npixels is the number of pixels in the current section of 
data being processed. The value of min_npixels is then used 
to calculate a merge acceleration factor, factor, which is mul- 
tiplied times the dissimilarity criterion value. If small_npix is 
the number of pixels in the smaller of the two regions being 
compared, factor 1 .0 if small_npix^min_npixels and 


factor=1.0-((min_npixels-small_npix)/min_np ixels), 

otherwise. 


35 


40 


(26b) 


45 


50 
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seam_threshold_factor (float) 


This threshold factor is used in 
determining whether a region found 
across a processing window seam is to 
be considered in determining whether a 
pixel is to be split out of its current 
region. 

(1 .0 <= seam_threshold_factor, 
default = 1.5. If threshold_factor = 

1.0, no regions are selected by this 
method) 


During the processing window elimination process, a “can- 
didate region label” set is accumulated for use in considering 
whether or not a pixel is to be split out of its current region. A 
candidate region is a region that either may contain pixels that 
should be split out and possibly be assigned to a different 
region, or is a region to which a split out pixel may be assigned 
to. Consider the data points that are in the pairs of rows and 
columns along the seam between the data quadrants reas- 
sembled in step 2 of the RHSEG algorithm. For each of these 
pixels calculate the dissimilarity between the pixel and its 
current region (own_region_dissim), and calculate the dis- 
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similarity between the pixel and the region of the pixel across 
the seam (other_region_dissim). If own_region_dissim> 
s eam_thresho ld_factor !: o ther_region_di s s im, add the region 
label of the region of the pixel across the seam to the “candi- 
date region label” set of the region the pixel belongs to. 


region_threshold_factor (float) This threshold factor is used in 

determining which regions are to be 
considered in determining whether a 
pixel is to be split out of its current 
region. 

(1.0 < threshold_factor, default =1.5. 
If region_threshold_factor = 

1.0, no regions are selected by 
this method) 


During the processing window elimination process, a “can- 
didate region label” set is accumulated for use in con- 
sidering whether or not a pixel is to be split out of its 
current region. Compare each region to every other 
region. If the dissimilarity between a pair of regions is 
less than region_threshold_factor*maxJhreshold, add 
each region label to the “candidate region label” set for 
the other region. NOTE: max_threshold is the maximum 
merging threshold encountered in the previous merging 
iterations. 


split pixels factor Pixel splitting factor. A pixel will be split 

(float) out from its current region if it is this 

factor more similar to another 
region than it is to its current region. 

(0.0 == split pixels factor, 

default = 1.5. No pixel splitting is 
performed if split pixels factor < 1 .0.) 


For each region with a non-empty “candidate region label” 
set, compute the dissimilarity of each pixel in that region 
to its current region (own_region_dissim) and to each 
region in the region’s “candidate region label” set (oth- 
er_region_dissim). If a pixel is found to have 
own_region_dissim>split_pixels_factor*other_region_ 
dissim, the pixel is split out from its current region. 


conv__nregions(short Number of regions for final convergence 

unsigned int) (0 < conv_nregions <65535, default = 2) 


When spclust_wght>0.0, the following optional parameters 
may be used to output information on closed connected 
regions (no defaults, ignored if spclust_wght=0.0): 


nb_conn_regons (string) Output number of connected regions list file 
conn_oparam (string) Connected regions output parameter file name 
(ignored if conn_rlblmap and conn_rnpixlist 
are not provided) 


Similar to the output parameter file oparam. 


conn_rlblmap (string) Output connected region label map data file name 
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Similar to rlblmap, but with closed connected regions. 


conn_regmerges (string) Output connected region label map merges list 
file name 


Similar to regmerges, but for closed connected regions. 


conn_mpixlist (string) Output connected region number of pixels list 
file name 


Similar to mpixlist, but for closed connected regions. 


conn_rmeanlist (string) Output connected region mean list file name 


Similar to nneanlist, but for closed connected regions. 


conn_rthreshlist (string) Output connected region maximum merge 
threshold list file name 


Similar to rthreshlist, but for closed connected regions. 


conn boundary_npix (string) Output connected region number of 

boundary pixels list file name 


Similar to boundary_npix, but for closed connected 
regions. 


conn_reg std dev (string) Output connected region standard 

deviation list file name 


Similar to reg_std_dev, but for closed connected regions. 

The following optional parameters control the run-time 
screen and log file outputs: 


debug (short unsigned int) Debug option 

(debug 0, default =1) 


Must be specified if debug>0 (ignored if debug=0): 


log_file (string) Output log file 


At a minimum (for debug=l), the output log file records 
program parameters and the number of regions and 
maximum merge ratio value for each level of the region 60 
segmentation hierarchy. 

The parameters that have the most effect on the nature of 
the segmentation results are spclust_wght, dissim_crit, 
chk_nregions and convfact. The default values are recom- 
mended for the other optional parameters for routine use of 65 
HSEG and RHSEG, with the exception that specification of 
the output file name parameters regmerges, rthreshlist, 


32 

boundary_npix and boundary_map is also recommended. In 
addition, the file name parameter rmeanlist is recommended 
when the number of spectral bands is less than 10 or so. Of 
course, if some input data elements are invalid, the some 
5 method of data masking should also be employed. 

The following paragraphs give some guidance on the set- 
ting of the spclust_wght, dissim_crit, chk_nregions and con- 
vfact parameters: 

10 spclust_wght: The user may want to vary the value of 
spclust_wght to modify the overall nature of the segmentation 
results. For spclust_wght=0.0, obtains relatively coherent 
closed connected regions. For spclust_wght=l .0, obtains 
5 relatively variated regions consisting of possibly several spa- 
tially disjoint subsections. For other values of spclust_wght 
obtains results intermediate the spclust_wght=0.0 and 
spclust_wght=1.0 results. 

dissim_crit: The user may also want to vary the value of 
20 dissim_crit to modify the overall nature of the segmentation 
results. The different dissimilarity criterion will result in dif- 
ferent merge ordering. NOTE: If any of the vector norm 
dissimilarity criterion is chosen (selections 1 , 2 or 3), the user 
25 may also want to specify the value of min_npixels_pct (small 
values in the range of 0. 1 to 1.0 are suggested). 

chk_nregions: The user may want to vary the value of 
chk_nregions to vary the level of segmentation detail in the 
most detailed level of the segmentation hierarchy. Higher 
30 values will increase the detail (the segmentation will have 
more regions) and lower values will decrease the detail (the 
segmentation will have fewer regions). 

convfact: The user may want to vary the value of convfact 
35 to control the number of hierarchical segmentation levels 
contain in the final segmentation hierarchy. Lower values of 
convfact will produce more segmentation levels and higher 
values of convfact will produce fewer segmentation levels. 
Too high of a value for convfact will produce only two hier- 
40 archical levels : one with chk_nregions number of regions and 
one with conv_nregions number of regions. 

Varying the other optional parameter values away from the 
default values requires a thorough understanding of the inner 
45 workings of the software implementation of the HSEG and 
RHSEG algorithms. 

This disclosure has been written in a general way as to 
include alternate embodiments of the algorithms in software. 
Varying the parameters will bring about various alternate 
50 embodiments of the innovation. 

The extension of the basic RHSEG algorithm to 3 -spatial 
dimensions is conceptually straightforward. One just needs to 
deal in data “voxels” instead of “pixels,” extend the definition 
of conn-type to include 3 -spatial dimension neighborhoods, 
and recursively divide the data into eight equal sections (halv- 
ing each spatial dimension) rather than four. 

Concerning conn_type: 


conn_type (short unsigned int) Neighbor connectivity type 
(3 -spatial dimensions) 

1. “Six Nearest Neighbors,” 

2. “Eighteen Nearest Neighbors,” 

3. “Twenty-Six Nearest Neighbors,” 
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For 2-spatial dimensions, we defined neighbors based on the 
following neighborhood chart, where the focal pixel is 
marked “X”: 
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Using this chart, N Nearest Neighbors include pixels 1, 

2, . . . N for 2-spatial dimension data. A similar chart could be 
developed for 3 -spatial dimensions, but printing would 15 
require 3 -dimensional plotting. It is easier to present a table of 
values. For 2-spatial dimensions, the table corresponding to 
above 2-spatial dimension chart, for up to Eight Nearest 
Neighbors, is: 


Neighbor 

Number 

nbcol 

nbrow 

1 

col - 1 

row 

2 

col + 1 

row 

3 

col 

row - 1 

4 

col 

row + 1 

5 

col - 1 

row - 1 

6 

col + 1 

row + 1 

7 

col + 1 

row - 1 

8 

col - 1 

row + 1 


where the focal pixel is at location (col, row), andnbcol is the 
neighboring pixel column location and nbrow is the neigh- 35 
boring pixel row location. 

The corresponding chart for 3 -spatial dimensions, for up to 
Twenty-Six Nearest Neighbors, is given on the next page. In 
this chart, the focal pixel is at location (col, row, depth), nbcol 
is the neighboring voxel column location, nbrow is the neigh- 40 
boring voxel row location and nbdepth is the neighboring 
voxel depth location. Note that the equivalent of Eight Near- 
est Neighbors (the default) in 2-spatial dimensions is Twenty- 
Six Nearest Neighbors in 3-spatial dimensions. 

The extension to 3-spatial dimensions of RHSEG is also 45 
conceptually straightforward. The seams between the data 
sections that are reassembled after processing at the previous 
level of recursion now become pairs of surfaces rather than 
pairs of lines. 


Neighbor 

Number nbcol nbrow nbdepth 


1 

col - 1 

2 

col + 1 

3 

col 

4 

col 

5 

col 

6 

col 

7 

col - 1 

8 

col + 1 

9 

col + 1 

10 

col - 1 

11 

col 

12 

col 

13 

col 

14 

col 

15 

col - 1 


row 
row 
row - 1 
row + 1 
row 


row 
row - 1 
row + 1 
row - 1 
row + 1 
row - 1 
row + 1 
row + 1 



depth 
depth 
depth 
depth 
depth - 1 
depth + 1 
depth 
depth 
depth 
depth 
depth - 1 
depth + 1 
depth - 1 
depth + 1 
depth - 1 
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-continued 


Neighbor 

Number 

nbcol 

nbrow 

nbdepth 

16 

col + 1 

row 

depth + 1 

17 

col + 1 

row 

depth - 1 

18 

col - 1 

row 

depth + 1 

19 

col - 1 

row - 1 

depth - 1 

20 

col + 1 

row + 1 

depth + 1 

21 

col + 1 

row + 1 

depth - 1 

22 

col + 1 

row - 1 

depth + 1 

23 

col + 1 

row - 1 

depth - 1 

24 

col - 1 

row + 1 

depth + 1 

25 

col - 1 

row + 1 

depth - 1 

26 

col - 1 

row - 1 

depth + 1 


The extension of the parallel implementation of RHSEG to 
3-spatial dimensions is also straightforward conceptually. 
The number of processors utilized by the 3-spatial dimension 
implementation is numprocs=8*”^ eve * y-1 . At inb_levels=2, 
numprocs=8, at inb_levels=3, numprocs=64, and at inb_lev- 
els=4, numprocs=5 1 2 . 

At program initialization, the input data would again be 
parceled out to numprocs processors, and each processor 
would independently process each of the numprocs data sec- 
tions. Again similar to the 2-spatial dimension version, upon 
completion of processing at recursive levels mb_levels 
through inb_levels, the 3-spatial dimension version would 
transfer its results to the appropriate processor at recursive 
level inb_levels-l. 

The parallel implementation scheme for the 3-spatial 
dimension version of RHSEG could be charted similarly to 
FIG. 14 , except that each recursive level splits eight ways 
rather than four ways. 

While the parallel implementation of the modified version 
of RHSEG is noted to be conceptually straightforward, the 
amount of detail that needs to be kept track of by the program 
is increased substantially, complicating the programming 
task substantially. 

All publications, patents, and patent documents are incor- 
porated by reference herein as though individually incorpo- 
rated by reference Although preferred embodiments of the 
present invention have been disclosed in detail herein, it will 
be understood that various substitutions and modifications 
may be made to the disclosed embodiment described herein 
without departing from the scope and spirit of the present 
invention as recited in the appended claims. 
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I claim: 

1. A computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of eliminating processing window artifacts that occur 
in segmentation of data having spatial characteristics through 
a recursive approximation of a segmentation process using a 
computer to perform the steps of: : (a) recursively dividing the 
data into subsections, each subsection having a boundary and 
no more than a predetermined number of data elements, des- 
ignating as contagious each new region that is adjacent to a 
data subsection boundary; (b) calculating a dissimilarity cri- 
terion between each new region and each other spatially adja- 
cent region; (c) if the dissimilarity criterion between one new 
region and another spatially adjacent region is less than or 
equal to a maximum merging threshold and if the one new 
region is not contagious, merge the one new region with a 
predetermined spatially adjacent region, further including if 
the predetermined number of merged regions in a data sub- 
section has not been reached, remove a contagious designa- 
tion from all regions in that data subsection and compute a 
new dissimilarity criterion between each new region and each 
other spatially adjacent region and each candidate region: (d) 
repeat step (c) until a predetermined number of merged 
regions is attained or until no further merges can occur, (e) 
compare each region to every other region through said dis- 
similarity criterion, and if said dissimilarity between a pair of 
regions is below a predetermined value, designate each of the 
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pair of regions a candidate region of the other and, further 
including for each pixel designated as a split pixel, associate 
candidate regions to the pixel equal to the candidate regions 
associated with the region originally assigned to the pixel, 
5 adding the original region assignment of the pixel to the 
candidate regions for the pixel, and the split pixel out of the 
region to which the pixel was originally assigned. 

2. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 

10 method of eliminating processing window artifacts according 
to claim 1, further including if the compared pair of regions 
has dissimilarity less than a spatial weighting factor less than 
the maximum merging threshold, merge the pair of regions. 

3. The computer implemented singleband multispectral 
15 remotely sensed imagery data for earth science observation 

method of eliminating processing window artifacts according 
to claim 2, further including returning a segmentation result 
to a calling recursive level program element. 

4. The computer implemented singleband multispectral 
20 remotely sensed imagery data for earth science observation 

method of eliminating processing window artifacts according 
to claim 1, further including, for each region with one or more 
candidate regions, computing a dissimilarity criterion for 
each pixel to the region and to each of the candidate regions. 
25 5. The computer implemented singleband multispectral 

remotely sensed imagery data for earth science observation 
method of eliminating processing window artifacts according 
to claim 4 , further including if the dissimilarity of the pixel to 
the region of the pixels across the boundary, when weighted 
30 by a predetermined factor, is less than the dissimilarity of the 
pixel to its own region, then designate each of the pair of 
regions as a candidate region of the other. 

6. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 

35 method of eliminating processing window artifacts according 
to claim 5 , further including for each region with one or more 
candidate regions, compute a dissimilarity criterion for each 
pixel to the region and to each of the candidate regions, and if 
the dissimilarity criterion for a pixel is less than the dissimi- 
40 larity criterion for one of the candidate regions when 
weighted by a predetermined factor, then designate said pixel 
as a split pixel. 

7. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 

45 method of eliminating processing window artifacts according 
to claim 1, wherein the dissimilarity criterion is selected from 
the group of vector norms, the square root of band sum mean 
squared error, the square root of band maximum squared 
error, and the SAR Speckle Noise Criterion. 

50 8. A computer implemented singleband multispectral 

remotely sensed imagery data for earth science observation 
method of eliminating processing window artifacts occurring 
in a segmentation of data with spatial characteristics compris- 
ing the steps of: (a) implementing a recursive segmentation 
55 algorithm on the data thereby assigning data pixels compris- 
ing the data into regions, in a fashion that divides the data into 
subsections; (b) for pixels along a seam between data subsec- 
tions determining whether a first pixel is within a predeter- 
mined criteria to its current region or to the region of a second 
60 pixel across the seam; (c) if the first pixel is determined to be 
within a predetermined criteria to the region of the second 
pixel across the seam by a predetermined threshold factor, 
then the second pixels region is designated as a candidate 
region for assignment of the first pixel; (d) comparing a first 
65 region to every other region within a predetermined threshold 
and if the similarity between a pair of regions exists then 
designating each of the pair of regions as a candidate region of 
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the other; (e) for each first region with a candidate region, 
computing the similarity of pixels in the first region to the first 
region and to candidate regions of the first region; (f) if a pixel 
in the first region is within a predetermined criteria to a 
candidate region, designating said pixel as a split pixel; (g) if 5 
merges between spatially adjacent and non-adjacent regions 
have equal importance then meting the pixels designated as 
split pixels into the regions to which those split pixels are 
within said predetermined criteria: (h) if merges between 
spatially adjacent regions have importance within said pre- 10 
determined criteria than non-adjacent regions then designat- 
ing each split pixel adjacent to a seam as contagious; (i) if a 
first split pixel not designated as contagious is within a pre- 
determined criteria to a candidate region not designated as 
contagious, merging the first split pixel into that most similar 15 
candidate region; (j) if either the first split pixel or the pixel 
within a predetermined criteria candidate region is designated 
as contagious, then designating both as contagious and defer- 
ring consideration of their merger at the present level of 
recursion; (k) returning to step (b) at a higher level of recur- 20 
sion until the top recursion level or a small number of regions 
is achieved. 

9. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of claim 8 wherein the recursive segmentation algo- 25 
rithm includes a region growing approach. 

10. The computer implemented singleband multispectral 

remotely sensed imagery data for earth science observation 
method of claim 8 wherein the recursive segmentation algo- 
rithm includes a special clustering approach. 30 

11. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of claim 8 wherein the data with spatial characteris- 
tics is image data. 

12. The computer implemented singleband multispectral 35 
remotely sensed imagery data for earth science observation 
method of claim 8 wherein the determination of the similarity 
of a pixel or region utilizes a dissimilarity criterion. 

13 . The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 40 
method of claim 8 wherein the dissimilarity criterion is 
selected from the group of vector norms, the square root of 
band sum mean squared error, the square root of band maxi- 
mum squared error, and the SAR Speckle Noise Criterion. 

14. A computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of eliminating window artifacts during the processing 
of a data set of pixels with spatial characteristics, wherein the 
data set is recursively divided into subsections and segmented 
into a plurality of regions through a recursive segmentation 
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algorithm using a computer to perform the steps of beginning 
below a top level of recursion (a) identifying first pixels along 
boundaries of the subsections at the current level of recursion; 
(b) designating the identified first pixels that are within a 
predetermined criteria to a second pixel along an adjacent 
boundary of another subsection as having a candidate region 
set including the second pixel’s region; (c) setting each iden- 
tified first pixel with a non-empty candidate region set as a 
separate region and designating each such region as a conta- 
gious region; (d) for each first region computing the dissimi- 
larity of the first region with every other region and if the 
dissimilarity is less than a threshold value then including such 
other similar region in a candidate region set for the first 
region, wherein in computing the dissimilarity of first regions 
with every other region, spatially adjacent regions are 
weighted to be similar within a predetermined criteria than 
equivalent non-adjacent regions: (e) if a first region with a 
non-empty candidate region set is a contagious region then 
designating within a predetermined criteria region of that first 
region’s candidate region set as contagious; (f) at a top recur- 
sion level, remove all contagious designations from regions; 
(g) if a first region is not contagious and the most similar 
region in its candidate region set is not contagious, then 
merging the first region and most similar candidate region; (h) 
unless at the top recursion level, returning to step (a) at a 
higher level of recursion until the number of regions is less 
than a predetermined number. 

15. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of claim 14 wherein the recursive segmentation algo- 
rithm includes a region growing approach. 

16. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of claim 14 wherein the recursive segmentation algo- 
rithm includes a special clustering approach. 

17. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of claim 14 wherein the data with spatial character- 
istics is image data. 

18. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of claim 14 wherein the determination of the similar- 
ity of a pixel or region utilizes a dissimilarity criterion. 

19. The computer implemented singleband multispectral 
remotely sensed imagery data for earth science observation 
method of claim 18 wherein the dissimilarity criterion is 
selected from the group of vector norms, the square root of 
band sum mean squared error, the square root of band maxi- 
mum squared error, and the SAR Speckle Noise Criterion. 



